% estimate impulse responses to government spending shock 
% forecast error identification
% this version: May 2011

clear all; clc; close all;

load us_quarterly       % generated by getdata.m


saveirf      = 1;

rand('state',1);
nptsvar     = 30;       % number of points of IRF in VAR
teststatio  = 0;        % tests wether VAR is explosive
confidence  = 1.645;    % (1.645: 90 percent);
ndraws      = 1000;      % number of replications in bootstrap
nlagq       = 4;        % number of lags
beg_samp    = (82-47)*4+1; 
end_samp    = (107-47)*4+4;

disp(['Sample: ' num2str(sampleq(beg_samp+nlagq)) ' - '  num2str(sampleq(end_samp)) ' (dependent variables)  '  ])
gy = mean(g(beg_samp:end_samp)./y(beg_samp:end_samp));
cy = mean(c(beg_samp:end_samp)./y(beg_samp:end_samp));
iy = mean(inv(beg_samp:end_samp)./y(beg_samp:end_samp));

    
for jrun = 1:3;
    
   
   if jrun==1;
       namshort = char('news','g','y','c','lr','rx','dy','inf');       
       xdata = [news log(g) log(y) log(c) lrr log(rx) dy inflation];        
       yshare = [1 gy 1 cy 1 1 1 1];
    elseif jrun ==2;
        namshort = char('news','g','y','i','lr','rx','dy','inf');       
        xdata = [news log(g) log(y) log(inv) lrr log(rx) dy inflation];        
        yshare = [1 gy 1 iy 1 1 1 1];
    elseif jrun ==3;
        namshort = char('news','g','y','nx','lr','rx','dy','inf');       
        xdata = [news log(g) log(y) nx  lrr log(rx) dy inflation];        
        yshare = [1 gy 1 1 1 1 1 1];
    end

    
    %% specify deterministics
    [obs nvar]=size(xdata);
    exo = []; 
    for ik=1:obs; ltrend(ik)=ik; end
    for ik=1:obs; exosq(ik)=ik^2; end
    exo = [ones(obs,1) ltrend' ];
    exo = exo(beg_samp:end_samp,:);
    samplelength=length(xdata(beg_samp:end_samp,:));

    %% estimate impulse responses
    [impzout,erzout,azeroout,a0betazout,var_explode]=...
        estimateIRF(xdata(beg_samp:end_samp,:),nlagq,nptsvar,teststatio,exo);       

    
    if jrun ==1;
        
        irf1 =  impzout(:,:,1);
        irf2 =  impzout(:,:,2);
        
        %% bootstrap IRF
        simdout=simdata(xdata(beg_samp:end_samp,:),...
        nlagq,nptsvar,ndraws,impzout,erzout,azeroout,a0betazout,exo);
        simimpzout=[];
        jjk=1;
        for zx=1:ndraws;
            
            
            [simimpzout,erzout,azeroout,a0betazout,var_explode] = estimateIRF(simdout(:,:,zx),...
            nlagq,nptsvar,teststatio,exo);                   
            irf1_mc(jjk,:,:)=  squeeze(simimpzout(:,:,1));
            irf2_mc(jjk,:,:)=  squeeze(simimpzout(:,:,2));
            jjk=jjk+1;
        
        end
        irfse1 = squeeze(std(irf1_mc(:,:,:),0,1));
        irfse2 = squeeze(std(irf2_mc(:,:,:),0,1));
        
    elseif jrun ==2;
        irf1 = impzout(:,:,1);
        irf2 = impzout(:,:,2);

        %% bootstrap IRF
        simdout=simdata(xdata(beg_samp:end_samp,:),...
        nlagq,nptsvar,ndraws,impzout,erzout,azeroout,a0betazout,exo);
        simimpzout=[];
        jjk=1;
        for zx=1:ndraws;
            
            
            [simimpzout,erzout,azeroout,a0betazout,var_explode] = estimateIRF(simdout(:,:,zx),...
            nlagq,nptsvar,teststatio,exo);                   
            irf_mc21(jjk,:,:) =  squeeze(simimpzout(:,:,1));
            irf_mc22(jjk,:,:) =  squeeze(simimpzout(:,:,1));
            jjk=jjk+1;
        
        end
        irfse1 = squeeze(std(irf_mc21(:,:,:),0,1));
        irfse2 = squeeze(std(irf_mc22(:,:,:),0,1));
        
   elseif jrun ==3;
        irf1 = impzout(:,:,1);
        irf2 = impzout(:,:,2);

        %% bootstrap IRF
        simdout=simdata(xdata(beg_samp:end_samp,:),...
        nlagq,nptsvar,ndraws,impzout,erzout,azeroout,a0betazout,exo);
        simimpzout=[];
        jjk=1;
        for zx=1:ndraws;
            
            
            [simimpzout,erzout,azeroout,a0betazout,var_explode] = estimateIRF(simdout(:,:,zx),...
            nlagq,nptsvar,teststatio,exo);                   
            irf_mc21(jjk,:,:) =  squeeze(simimpzout(:,:,1));
            irf_mc22(jjk,:,:) =  squeeze(simimpzout(:,:,1));
            jjk=jjk+1;
        
        end
        irfse1 = squeeze(std(irf_mc21(:,:,:),0,1));
        irfse2 = squeeze(std(irf_mc22(:,:,:),0,1));
        
    end
    
    
    
    
    
    %% plot
    time = [0:nptsvar-1];
    jj=1;
    figure
    for zx=1:nvar;
        subplot(3,3,jj)
        mycolors = [1 1 1; .85 .85 .85];            
        colormap(mycolors);
        rescale = max(irf1(1:10,2))*(yshare(2)/yshare(zx));
        point = irf1(:,zx)./rescale;
        sepoint = irfse1(:,zx)./rescale;
        ub = point+confidence*sepoint;
        lb = point-confidence*sepoint;  
        hold on;   
        gx=[lb ub-lb];
        area(time,gx,min(lb),'EdgeColor','none');
        plot(time,zeros(nptsvar,1),'k--','LineWidth',1);
        plot(time,point,'b-','LineWidth',1);
        title(deblank(namshort(zx,:)));
        jj=jj+1;
    end
    
        if saveirf
        
        disp('writing figures to files...')
        for zx=1:nvar;            
            figure
            set(gcf,'Visible','off')
            mycolors = [1 1 1; .85 .85 .85];            
            colormap(mycolors);

            rescale = max(irf1(1:10,2))*(yshare(2)/yshare(zx));
            point = irf1(:,zx)./rescale;
            sepoint = irfse1(:,zx)./rescale;
            ub = point+confidence*sepoint;
            lb = point-confidence*sepoint;  

            hold on;   
            gx=[lb ub-lb];
            area(time,gx,min(lb),'EdgeColor','none');
            plot(time,zeros(nptsvar,1),'k--','LineWidth',1);
            plot(time,point,'b-','LineWidth',4);
            axis([-.2 nptsvar-.8  min([1.5*min(lb) 0]) 1.5*max(ub)]);box on;       
            set(gca,'FontSize',40)
            box on;
            saveas(gcf,['..\figures\Ram_' num2str(jrun) '_' deblank(namshort(zx,:)) '.eps'],'psc2')            
        end

    end

end